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Abstract 

Axions are one of the most attractive dark matter candidates. The evolution of their number density in the early uni¬ 
verse can be determined by calculating the topological susceptibility xiT) of QCD as a function of the temperature. 
Lattice QCD provides an ab initio technique to carry out such a calculation. A full result needs two ingredients; physi¬ 
cal quark masses and a controlled continuum extrapolation from non-vanishing to zero lattice spacings. We determine 
X{T) in the quenched framework (infinitely large quark masses) and extrapolate its values to the continuum limit. 
The results are compared with the prediction of the dilute instanton gas approximation (DIGA). A nice agreement 
is found for the temperature dependence, whereas the overall normalization of the DIGA result still differs from the 
non-peiturbative continuum extrapolated lattice results by a factor of order ten. We discuss the consequences of our 
findings for the prediction of the amount of axion dark matter. 
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1. Introduction 


One of the greatest puzzles in particle physics is the 
nature of dark matter. A prominent particle candidate 
for the latter is the axion A [1, 2]: a pseudo Nambu- 
Goldstone boson arising from the breaking of a hypo¬ 
thetical global chiral t/(l) extension [3] of the Standard 
Model at an energy scale much larger than the elec- 
troweak scale. 

A key input for the prediction of the amount of axion 
dark matter [4—6] is its potential as a function of the 
temperature, V{A, T). It is related to the free energy 
density in QCD, F(0, T) = - In Z(6, T)I'V, via 


V(A,T) = --\n 


Z(0, T) 
Z(0, T) 


lfl=A/X,, 


( 1 ) 


where 'V is the Euclidean space-time volume. The an¬ 
gle 0 enters the Euclidean QCD Lagrangian via the ad¬ 
ditional term involving the topological charge density 


q{x), 

-i0q{x) = -i0^ef,ypc-F“,^{x)F‘‘pJx), (2) 

with F°y being the gluonic field strength and s 
/(Ttt) the fine structure constant of strong interactions. 
On general grounds, the free energy density and thus 
the axion potential has an absolute minimum at 0 - 
AI/a - 0. In fact, this is the reason why in this exten¬ 
sion of the Standard Model there is no strong CP prob¬ 
lem [3]. The curvature around this minimum determines 
the axion mass niA at finite temperature. 


mliT) = 


dMA,T)^ X{T) 

QA2 1^=0 fl 


(3) 


in terms of the topological susceptibility, i.e. the vari¬ 
ance of the 0 — 0 topological charge distribution, 

X{T) = f d'^x(q(x)q(0))Tls^o = Hm ( 4 ) 

I ' k —>00 y 
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where Q = J d^xq(x). Similarly, self-interaction terms 
in the potential, e.g. the term occurring in the expan¬ 
sion of the free energy density around 6 = A//a - 0, 


y(A, T) = [l + b 2 (T)(^ + ...] |fl=A//„ (5) 


are determined by higher moments. 


biiT) = 


- 3<e2)2 


( 6 ) 


These non-perturbative quantities enter in a prediction 
of the lower bound on the fractional contribution of ax- 
ions to the observed cold dark matter as follows* 


Ra > f({e^)) 


GeV^ 

g*s(Tosc)Tasc 



(7) 


where {6^) is the variance of the spatial distribution 
of the initial values of the axion field Ajf a before the 
formation of the dark matter condensate by the mis¬ 
alignment mechanism, which occurs, when the Hub¬ 
ble expansion rate gets of the order of the axion mass, 
ifiAiTosc) - 3//(rosc), i-e. at a temperature 


50 GeV 


( 8 ) 


The function f{{6^}) is taking into account anharmonic- 
ity effects arising from the self-interaction terms in 
V(A, T) and depends on the specific form of the poten¬ 
tial. The functions gtR{T) and gts(T) denote the effec¬ 
tive number of relativistic energy and entropy degrees 
of freedom, respectively. 

What is urgently needed for axion cosmology is thus 
a precise determination of the topological susceptibility 
and higher moments of the topological charge distribu¬ 
tion. In this context, most predictions have been entirely 
based on the semi-classical expansion of the Euclidean 
path integral of finite temperature QCD around a dilute 
gas of instantons - finite action minima of the Euclidean 
action with unit topological charge - see e.g. Ref. [8] 
for an early exhausting study and Ref. [9] for a recent 
update concerning the quark masses. A comparative 
study of these predictions based on the dilute instanton 


* This is a rewriting of equation (2.10) of Ref. [7] where we have 
used;r(0) = (”M/a) 2 - 3.6 x 10“^ GeV"* from the chiral Lagrangian 
to express fA in terms of ma, the zero temperature mass, that is itself 
a function of T^sc and the ratio^(7’)/^(0) = (mA(T)! ma)^ through the 
condition for the onset of the oscillations, mAiTasc) = ^H(Tasc) with 
= 87T^GNgtR(T)T^/90 the Hubble expansion rate. 


gas approximation (DIGA) has been carried out in Ref. 
[10], where also an analysis in terms of a phenomeno¬ 
logical instanton liquid model (IILM) [11] is presented. 
However, up to now, in all DIGA investigations of the 
topological susceptibility only the one-loop expression 
in the expansion around the instanton background field 
was used. This results in a strong renormalization scale 
dependence and thus large uncertainties which were ne¬ 
glected in the previous DIGA based predictions. In fact, 
in the temperature range T^sc ~ GeV of interest, one 
expects a large uncertainty in the overall normalization 
due to the neglection of higher order loop effects, since 
in this region asiTosc) is not small. We will exploit 
in this letter both the one-loop DIGA result as well as 
its two-loop renormalization group improved (RGI) ver¬ 
sion in order to study the theoretical uncertainties aris¬ 
ing from higher loop corrections. Most importantly, we 
compare these predictions with the outcome of our lat¬ 
tice based fully non-perturbative results. 

Actually, there have been a number of lattice cal¬ 
culations of xiT) and b 2 (T) at temperatures below or 
slightly above the QCD phase transition, mostly in 
quenched QCD, see e.g. [12-14]. Here we go beyond 
those lattice calculations. We extend the available tem¬ 
perature range and carry out a controlled continuum ex¬ 
trapolation for this extended range. In addition, note 
that there has been no quantitative investigation whether 
and where the lattice results turn into the DIGA results. 
We will present our new high-quality lattice data for the 
topological susceptibility in quenched QCD (i.e. ne¬ 
glecting the effects of light quarks) and compare them 
quantitatively to the DIGA result. 

2. Axion potential coefficients from the lattice 

On the lattice, the topological susceptibility is mea¬ 
sured on the torus as the second moment of the distribu¬ 
tion of the global topological charge 

A/ = 

where Q is any renormalized discretization of the global 
topological charge, and A' is the four-volume of the lat¬ 
tice. There are a lot of different fermionic and gluonic 
definitions of Q available. We choose a gluonic defini¬ 
tion based on the Wilson flow [15, 16], which has the 
correct continuum limit similarly to the fermionic defi¬ 
nitions but is numerically a lot cheaper. In particular we 
evolve our gauge field configurations with the Wilson 
plaquette action to a flow time t and define the global 
topological charge as the integral over the clover defi¬ 
nition of the topological charge density. This definition 
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Figure 1: Demonstration that volume is sufficiently large to have neg¬ 
ligible finite volume corrections on The data for TjTc = 2 is 
multiplied by 3 for better visibility of the comparison. 


gives a properly renormalized observable when the flow 
time t is fixed in physical units. 

We use a tree-level Symanzik improved gauge ac¬ 
tion. Our temperatures range from below Tc up to ATc- 
Here is the critical temperature, which is the quan¬ 
tity used for scale setting. The critical temperatures 
for different lattice spacings were determined in ear¬ 
lier work [17, 18]. For the whole temperature range we 
keep the spatial lattice size approximately at L = 2/7’^. 
We checked with dedicated high volume runs at l.STc 
and 2Tc that this volume is sufficiently large to have 
negligible finite volume corrections on Q^. This is 
shown in Fig.l. Our spatial geometry is L x L x 2L 
to enable tests of subvolume methods which will be 
reported separately - here we only use the full vol¬ 
ume. For all temperatures we have three lattice spac¬ 
ings (aT) ' = 5,6,8 to be able to perform an inde¬ 
pendent continuum extrapolation for every temperature. 
The local heatbath/overrelaxation algorithm is used for 
the update, one sweep consists of 1 heatbath and 4 over¬ 
relaxation steps. We found that the autocorrelation time 
of the topological charge depends weakly on (aT), i.e. 
if the temperature is increased by decreasing the lattice 
spacing. The number of update sweeps between mea¬ 
surements was chosen in accordance with the autocor¬ 
relation time. Tab. 1 lists the simulation points with the 
number of sweeps. 

We integrated the Wilson flow numerically to a max¬ 
imum flow-time of about 8f 1/(27’^) for all tempera¬ 
tures. Figure 2 gives the dependence of the susceptibil¬ 
ity on the flow time for 7’ - 2Tc. While in the contin¬ 
uum limit the result is independent of the choice of the 
flow time f, different choices have very different lattice 
artefacts. For small flow times the different lattice spac- 


T/T, 

Nr 

Ns 

^sweeps 

0.9 

5 

12 

32K 


6 

12 

48K 


8 

16 

170K 

1.0 

5 

12 

48K 


6 

12 

64K 


8 

16 

180K 

1.1 

5 

12 

48K 


6 

16 

160K 


8 

20 

330K 

1.3 

5 

16 

64K 


6 

16 

220K 


8 

24 

550K 

1.5 

5 

16 

96K 


6 

20 

210K 


8 

24 

660K 

1.7 

5 

20 

260K 


6 

20 

300K 


8 

28 

700K 

2.0 

5 

20 

280K 


6 

24 

510K 


8 

32 

840K 

2.3 

5 

24 

420K 


6 

28 

530K 


8 

36 

900K 

2.6 

5 

28 

710K 


6 

32 

2600K 


8 

44 

5000K 

3.0 

5 

32 

810K 


6 

36 

7400K 


8 

48 

1700K 

3.5 

5 

36 

910K 


6 

44 

2100K 


8 

56 

5100K 

4.0 

5 

40 

810K 


6 

48 

2200K 


8 

64 

4900K 


Table 1: List of simulation points: temperatures, lattice sizes Nt 
\l{aT), Ns = Lja, and number of sweeps are given. 
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Figure 2: Demonstration that Wilson flow/lattice renormalization are 
under control. 


Figure 3: Lattice data on the topological susceptibility at N[ = 5,6, 8 
and lattice continuum extrapolation together with ht of simple power 
law 


ings give very different results. For larger flow times 
the expected plateau behavior can be observed for each 
lattice spacing and the lattice artefacts also decrease sig- 
nihcantly. The choice of the flow time brings in some 
arbitrariness into the analysis, however the continuum 
result should not depend on this choice once t is fixed in 
physical units. But this is certainly a subleading source 
of error compared to the statistical error due to the rare 
topology tunneling events at high temperatures. In this 
analysis we choose a temperature dependent flow time 
for the evaluation of as 


T<1.5 

jl/r^, T>\.5' 


(9) 


For low/high temperatures this means a temperature in¬ 
dependent/dependent flow time. This choice is safely in 
the expected plateau region for all temperatures. 

The resulting values for the susceptibility are plotted 
in Fig. 3. This plot also gives the result of a global 
continuum extrapolation using a set of temperatures and 
the 6-parameter power law ansatz 

( y, \b+h'a^ 

where Tq, and b are fit parameters giving the contin¬ 
uum limit, Tq, and b' are ht parameters describing 
the deviation from the continuum limit. The ht param¬ 
eter To is included as a consistency check and should 
give 1 in units of Tc- This is satished by the ht result. 
The variation between different choices for the starting 
temperature of the ht range T^inlTc - 1.3,1.5,1.7 gives 
an estimate of the systematic error of the result. The 


best ht parameters are 

^0=0.11(2X1), h =-7.1(4)(2), To^ 1.02(5)(2), 

( 11 ) 

where the hrst error is the statistical, the second is the 
systematic. The point-wise continuum extrapolation is 
consistent with the global ht, evidently the latter has 
smaller errors for large temperatures. Note that though 
a controlled continuum extrapolation is possible using 
three lattice spacings, estimating the systematic uncer¬ 
tainty of this extrapolation would require at least one 
more lattice resolution. 

In a recent analysis [13] the topological susceptibil¬ 
ity was calculated using the techniques [19, 20]. The 
calculation was carried out at two temporal extensions, 
corresponding to two lattice spacings at each temper¬ 
ature. The exponent b - -5.64(4) was found which 
differs from our value. Note however, that our temper¬ 
ature range is larger, thus we are closer to the applica¬ 
bility range of the DIGA. Furthermore, the two lattice 
spacings were not sufficient for a controlled continuum 
extrapolation thus uncertainties related to this hnal step 
are not included in the result of [13]. 

We have also determined the second important coeffi¬ 
cient b 2 of the axion potential, characterizing its anhar- 
monicity, by measuring the observable 

, <e"> - 


The result is plotted in Fig. 4. 
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Figure 4: Lattice data on the anharmonicity coefficient ^2 of tho axion 
potential compared to its DIGA prediction. The data points are shifted 
a bit horizontally for better visibility. 

3 . Comparison between lattice and DIGA results 


Figure 5: Prediction of the topological susceptibility in the DIGA: 
comparison between one-loop and two-loop RGI results. We used the 
four-loop expression for the running coupling in the modified mini¬ 
mal subtraction scheme as given in the appendix of Ref. [29] and the 

central value of Tc/A— = 1.26(7) as determined from the lattice 
'ms ' 

in Ref. [18]. 


In this section, we confront the lattice results with the 
ones obtained from the DIGA framework. For the sake 
of completeness let us collect first the available formulas 
for the latter [21-28]. At very high temperatures, far 
above the QCD phase transition, it makes sense to infer 
the 9 dependence of QCD from the DIGA, in which the 
partition function can be written as [22], 

Z{9, T) - V —L^z;'^"'-(7’)exp [Wim - «;)] , (12) 
«;■«/ ‘ 

where Z/ = Zj is the contribution arising from the ex¬ 
pansion of the path integral around a single instanton I 
(anti-instanton /). It follows directly that the potential 
has the form 

V{A,T)^X(T){l-cos0)\e^Aih, ( 13 ) 

from which one infers 

(14) 

This can be confronted right-away with our lattice re¬ 
sults, cf. Fig. 4. Similar to Ref. [12] we find that the 
prediction from the DIGA for b 2 is reached already at 
surprisingly low values of T/Tc > 1. 

The whole temperature dependence of the axion po¬ 
tential arises in the DIGA through the topological sus¬ 
ceptibility, which in this case is explicitly given by 

A(r)-- — -=2J dpD{p)G(npT), (15) 

in terms of the instanton size distribution at zero tem¬ 
perature, Dip), and a factor G{npT) taking into account 


finite temperature effects. The former is known in the 
framework of the semiclassical expansion around the 
instanton for small as(pr) ln(p Pr) and p mi(pr), where 
as is the strong coupling, pr is the renormalization scale 
and mi(pr) are the running quark masses. To one-loop 
accuracy, it is given by^ 

2n I 2n 

- exp- 




y.{pprf° [l + Q(aj;^(juQ)], 

with 

= ^e-«34122. (17) 

At finite temperature, electric Debye screening pro¬ 
hibits the existence of large-scale coherent fields in the 
plasma, leading to the factor [24, 25], 

G{x) = exp [-2)? - 18A(x)}, (18) 

with 

Mx) In [l -H {npTfl3\ -b a [l -b y{npTy^l^\ ^ , 

( 19 ) 

and a = 0.01289764 and y = 0.15858, in Eq. (15). This 
factor cuts off the integration over the size distribution 
in Eq. (15) at x = npT 1 and ensures the validity of 
the DIGA at large temperatures, at which asinT) ^ 1. 


^For quenched QCD the number of light quarks is tif = 0; the 
general formula can be found explicitly in e.g. Ref. [28], which con¬ 
tains a pioneering confrontation of cooled lattice data on D{p) with 
the two-loop RG improved DIGA result. 
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TITr 

1.5 

2 

3 

4 

5 

II 

o 

-6.04 

-6.26 

-6.43 

-6.50 

-6.55 

b(K^l) 

-6.37 

-6.46 

-6.55 

-6.59 

-6.62 

II 

-6.55 

-6.59 

-6.64 

-6.67 

-6.69 


Table 2: Temperature slopes of the topological susceptibility predicted 
in the two-loop RGI DIGA, for a range of renormalization scales ac¬ 
cording to Eq. (24). 


T/Tr 

1 

2 

3 

4 

5 


0.36 

0.23 

0.19 

0.17 

0.16 


1 

0.1 

0.01 

0.001 

0.0001 

1e-05 

1e-06 




1-loop |i=0 

1- loop |I=2 

2- loop |i=0 
2-loop |i=2 
iattioe oonti 

6/Pm . 

0 Pm 

6/d.^ 


O/Pm . 

luum I— 9- 








: 





ji,. 















Table 3: Strong coupling constant at fir = 11 pm for the temperature 
range covered by the lattice. 


Collecting all the factors, the topological susceptibil¬ 
ity, in the one-loop DIGA, reads 

\6 


X(T) 


2d^(nT)U^ 


I 

/ 2 ;r y 

[tttI 

1 ®Ms(^''^ 1 


( 20 ) 


X exp 


In 






with 


Jo 


dxx^G(x) ^0.267271. 


( 21 ) 


This result, however, still suffers from a sizeable de¬ 
pendence on the renormalization scale reflecting the 
importance of the neglected two-loop and higher order 
contributions. In fact, it is tamed by taking into account 
the ultraviolet part of the two-loop correction. The lat¬ 
ter has been calculated in Ref. [27] and shown to have 
exactly the form that the gauge coupling becomes a pa¬ 
rameter running according to the renormalization group 
(RG). Therefore, the ultimate, all order result for the 
topological susceptibility becomes independent of 
for jUr —> oo. At two loop, the corrections amount to 
a factor 


(p^^)0Si-12/?o)«HsO^A(4;r). ^ 102, (22) 


in Dip). Therefore, this RG improvement can be taken 
into account by replacing the factor I in Eq. (20) by 


I 


^ l^r j-30aHsO/,)/(4;r) 





G(x). 


(23) 


In fact, exploiting the two-loop RG improvement, the 
fir dependence is heavily reduced, as is obvious from 
Fig. 5, where we have used as a natural renormalization 
scale 

fir ^ kIp„^ kttTI 1.2, (24) 


Figure 6: Rescaled one-loop and two-loop RGI DIGA results com¬ 
pared to lattice continuum extrapolation. The DIGA results shown in 
Fig. 5 were scaled by a factor K of order ten such that they coincide 
at T jTc = 2 with the central value of the lattice continuum data. 


with Pm being approximately the maximum of the in¬ 
tegrand of 7, and varied the remaining free parameter 
K between 0.6 and 2. The renormalization scale depen¬ 
dence appears to be highly reduced in the regimes > 2>Tc 
and < Tc- However, in an intermediate region, ~ Tc- 
2Tc, it is comparable in size to the one at one-loop. 

We present in Table 2 the power-law behavior pre¬ 
dicted by the two-loop RGI DIGA at various tempera¬ 
tures, which can be compared to the fit (11) to the con¬ 
tinuum lattice result. As far as the overall normalization 
of the DIGA result for;^^ is concerned, one still expects 
a large uncertainty in the temperature range available 
from the lattice. In fact, at these temperatures, is not 
small, see Tab. 3. Apart from the ultraviolet part, there 
will be a finite part of the two-loop correction which will 
affect mainly the overall normalization of x and will de¬ 
pend on the temperature only logarithmically. Unfor¬ 
tunately, this finite part is not known, yet. Therefore, 
when comparing to the continuum extrapolated lattice 
results, we allow a multiplicative factor K to account 
for this uncertainty, i.e. we absorb the remaining higher 
loop uncertainties by replacing 

[l + 0(a^(fir = KnT/1.2))] ^ K(T/Tr) (25) 

in Eq. (20) and the corresponding two-loop RGI expres¬ 
sion. Clearly, the /T-factor should approach unity at very 
large r/Tc. 

Figure 6 nicely illustrates the agreement between the 
DIGA and the lattice result, if a /T-factor of order ten is 
included^. More precisely, fitting the lattice continuum 


factors of order ten to fifty are not uncommon even at next-to- 
leading order in ordinary perturbative QCD, see e.g. Ref. [30]. 
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data with the rescaled DIGA expression in the tempera¬ 
ture range TjTc > 1, one finds 

K^9.2± 0.6, at 95% CL, (26) 

while a fit in the temperature range T/Tc >2 yields 

8.0 ±3.5, at 95% CL. (27) 

Apparently, in the temperature range accessible to the 
lattice, the higher order corrections to the pre-factor of 
the DIGA are still appreciable, but there are indications 
of a trend that the ^f-factor gets smaller, as expected, 
towards larger values of T jTc. 

The /T-factor strongly depends on the value of 
Tcl it reduces to one for TdlA^—^^ ^ 1.03. 
However, the latter value is about 3 sigma below the 
central value determined in Ref. [18], TdlA—^'' ^ 

L J, Cl I 

1.26(7). 

4. Conclusions 

This paper presents lattice and DIGA calculations of 
the 6 and temperature dependence of the free energy 
density of QCD in the quenched limit, i.e. for infinite 
quark masses. In the lattice approach the temperature 
ranges from Q.9Tc to ATc and thus significantly extends 
former results. The precise high temperature data al¬ 
lows to compare the lattice results with first principle 
calculations within the DIGA. For the evaluation of 
X(T) we use the two-loop RGI version of the instanton 
size distribution, which is less sensitive with respect to 
the renormalization scale. We find that the two differ¬ 
ent methods show nice agreement after one includes an 
overall correction factor in the DIGA results. This order 
ten correction is supposed to take into account higher 
loop contributions in the normalization, similar to so- 
called /T-factors in perturbative QCD processes. These 
are significant, since in the considered temperature re¬ 
gion the strong coupling constant is still large. Never¬ 
theless, the achieved results motivate an investigation of 
more realistic models involving light quarks. 

The possible consequences of our findings on the pre¬ 
dictions of the amount of axion dark matter can be read 
off from Fig. 7. First, we have plotted the regions where 
axion DM would be overproduced with respect to ob¬ 
servations and are thus excluded. The darker yellow 
region is excluded just by the axions produced from 
the misalignment mechanism, cf. Eq. (7), assuming a 
flat distribution of initial misalignment angles in the ob- 


Quenched QCD 

2-loop RGI DIGA T,^ 294 MeV. K= 8±3.5 (k= 0.6-2) (gray) 



10-1 1 10 
fosclGeV] 


Unquenched QCD 

2-loop RGI DIGA K = 1 (blue), AT = 9.22±0.6 (gray) {k= 0.6-2) 
IILM from Ref. [10] (dashed red) 



Figure 7: Consequence of our findings for axion dai'k matter from 
quenched (top) and full QCD (bottom). The dark yellow region is 
excluded because Ra > I even without a contribution from strings. 
The light yellow region indicates Ra > I when string contributions 
are included. In the dark green region axions even including strings 
give a small contribution (Ra < 0.1) to dark matter. The light green 
region indicates Ra <0.1 just from the misalignment mechanism. In 
order to compare the quenched results with the axion dark matter sce¬ 
nario we had to transform dimensionless quantities into dimensionful 
ones. This is not unambiguous since we do not have pure gauge the¬ 
ory in nature. Using different quantities to set the scale one gets 280 
-310 MeV for Tc. For illustrative puiposes we use Tc = 294 MeV 
which lies in the middle of the range and is a factor two larger than the 
full QCD transition temperature. Our lattice results are shown by the 
blue points and the two-loop RGI DIGA prediction by the gray band 
in the quenched Figure. In the unquenched case the blue and gray 
bands correspond to the two-loop RGI DIGA predictions with ^ = 1 
and K = 9.22 ± 0.6, respectively. The dashed red line shows the IILM 
prediction of Refs. [10, 11]. 
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servable universe^ 0 € [-7r,7r]. For the DIGA poten¬ 
tial, xiT){l - cos(0)), we calculate {6^}f({0^)) - 6.6 
from numerically evolving the axion field with differ¬ 
ent initial conditions. The exclusion extends to the light 
yellow region if we consider the DM axions produced 
from the decay of unstable axionic strings and domain 
walls according to Ref. [32]. The excluded region has 
to be compared with our lattice results (blue dots) and 
the two-loop DIGA calculation with the fitted /T-factor, 
for which we have set the scale = 294 MeV (twice 
the dynamical value) for illustration purposes, see Fig. 
7 (top) . Despite the arbitrariness we can derive a num¬ 
ber of interesting conclusions. First note that our results 
correspond to a region where axion DM does not of¬ 
fer a viable cosmology because > 1. However, our 
lattice results are too short in T/Tc only by a factor of 
2-3. Further developments from the lattice side could in 
principle lead to a direct estimation of the axion mass. 
Meanwhile, the DIGA result allows a controlled extrap¬ 
olation of our results over the interesting region, that 
can be used to set a constraint (here merely illustrative): 
25 fieV < mA < 525 peV assuming that at least ten 
percent of dark matter can be attributed to axions. 

For the final calculation of the QCD axion abun¬ 
dance as dark matter candidate one needs to include 
the light quarks. This complicates the lattice evalua¬ 
tion by far. However, the DIGA approximation can 
still be carried out and self-consistently takes into ac¬ 
count the running of the quark masses as well as their 
influence on the strong coupling constant. We present 
this preliminary result as a blue band in Fig. 7 (bot¬ 
tom). From what we learned in the quenched case, this 
DIGA result still misses an 0(10) constant factor but 
the lack of lattice data in this case prevents from fit¬ 
ting the expected Tt'-factor explicitly. We can tenta¬ 
tively estimate it to be the required factor that makes 
X(T)lx(0) = 1 at Tc- Since the transition in full QCD is 
a crossover [33], Tc is not unambiguous. Using the tran¬ 
sition temperature defined by the chiral susceptibility, 
Tc = 147(2)(3) MeV [34-36], we obtain K = 9.22 + 0.6. 
Very interestingly, this agrees with the IILM model cal¬ 
culation of Ref. [10, 11] (red-dashed line) at the temper¬ 
atures of interest. This factor of 9.22 and the previous 8 
of the quenched limits do not affect very strongly the ex¬ 
traction of values of the axion mass, owing to the strong 
T -dependence of xiT). Using K - 9.22 and again as¬ 
suming at least ten percent axion contribution to dark 
matter we get the range 40 jueV < mA < 930 peV 


the scenario where inflation homogenizes one particular value 
of 6 across our universe, a prediction of the axion DM abundance is 
not possible, except perhaps in an anthropic sense [4, 31] 


while using TiT = 1 we would get 50 peV < mA 
1100 /teV, i.e. only a ~ 20% correction for an (9(10) K 
factor uncertainty. 
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